'''
Created on Jun 21, 2011

@author: oabalbin

This script parses the cosmic database 
Reports a list of mutations by genes.
If a gene list is provided the report list is only for those genes.


Gene name       Accession Number        HGNC ID Sample name     ID_sample       ID_tumour       Primary site    Site subtype    Primary histology       Hist
ology subtype   Mutation ID     Mutation CDS    Mutation AA     Mutation Description    Mutation zygosity       Mutation NCBI36 genome position Mutation NCB
I36 strand      Mutation GRCh37 genome position Mutation GRCh37 strand  Pubmed_PMID
A1CF    ENST00000374001 24086   Br27P   1312984 1223893 central_nervous_system  brain   glioma  astrocytoma_Grade_IV    40879   c.976G>A        p.E326K Subs
titution - Missense     het     10:52245937-52245937    -       10:52575931-52575931    -       18772396
'''

import sys
import numpy as np
from collections import defaultdict,deque

def list_of_names(inputfile, pfield):
    '''
    '''
    genlist=deque()
    for line in inputfile:
        if line.startswith('#'):
            continue 
        fields = line.strip('\n').split('\t')
        genlist.append(fields[pfield])
    
    return genlist


def read_cosmic_mutation_file(ifile,outfile,genelist=False):
    '''
    Parses cosmic mutation file.
    0 Gene name
    1 Accession Number
    2 HGNC ID
    3 Sample name
    4 ID_sample
    5 ID_tumour
    6 Primary site
    7 Site subtype
    8 Primary histology
    9 Histology subtype
    10 Mutation ID
    11 Mutation CDS
    12 Mutation AA
    13 Mutation Description
    14 Mutation zygosity
    15 Mutation NCBI36 genome position
    16 Mutation NCBI36 strand
    17 Mutation GRCh37 genome position
    18 Mutation GRCh37 strand
    19 Pubmed_PMID
    
    Output a bed-file. positions are represented as -1 pos +1
    '''
    ifile=open(ifile)
    
    ofile=open(outfile,'w')
    pfile=open(outfile+'.allcosmic.pos','w')
    sfile=open(outfile+'.ActionableGeneCosmic.pos','w')
    fields=defaultdict()
    header = ifile.next().strip('\n').split('\t')
    for i,h in enumerate(header):
        fields[i]=h
    cols=[17,0,10,11,12,13,14,18]
    # Type of cosmic mutations not included for
    # homozygous determination
    mut_not_included=['Substitution - coding silent','Deletion - Frameshift',
                      'Insertion - Frameshift','Insertion - In frame',
                      'Deletion - In frame','Unknown']
    #names=[location, geneName, mutCDS, mutAA, mutDes, mutZyg, strand]    
    pos_already_seen=set()
    i=0
    for l in ifile:
        f = l.strip('\n').split('\t')
        location, geneName, mutID,mutCDS, mutAA, mutDes, mutZyg, strand = \
        f[cols[0]],f[cols[1]],f[cols[2]],f[cols[3]],f[cols[4]],f[cols[5]],f[cols[6]],f[cols[7]]
        
        if len(location.split(':'))<2:
            continue
        
        chr,pos=location.split(':')[0],location.split(':')[1]        
        start, end = pos.split('-')[0],pos.split('-')[1]
        print start, end
        if mutDes in mut_not_included or start!=end:
            continue

        
        info=['geneName=',geneName,';mutID=',mutID,';mutCDS=',mutCDS, ';mutAA=',mutAA, 
              ';mutDes=',mutDes, ';mutZyg=',mutZyg]
        info = ",".join(info).replace(',','')
        line = [str(chr),int(start)-1,int(end)]+[info]+['0']+[strand] #'chr'+
        pline=[chr,'@',start] #'chr',
        this_pos=",".join(pline).replace(',','')
        if this_pos not in pos_already_seen:
            pos_already_seen.add(this_pos)
            
            if len(mutCDS.split('>'))==2:
                base_reference,base_alternative=mutCDS.split('>')[0][-1],mutCDS.split('>')[1]
                #pline=['chr',chr,'@',start,'|%s>%s'%(base_reference,base_alternative)]
                #pline=['chr',chr,'@',start]
                pfile.write(",".join(pline).replace(',','')+',')
            #print ",".join(line).replace(',','\t')
            if genelist and (geneName in genelist):
                ofile.write(",".join(map(str,line)).replace(',','\t')+'\n')
                sfile.write(",".join(map(str,[line[0],line[1]])).replace(',','\t')+'\n')
        else:
            continue
        
    
    ofile.close()
            

pfield=0
ifile='/data/databases/cosmic/v56/CosmicMutantExport_v56_151111.tsv'
ofile='/exds/users/oabalbin/projects/exomes/genelists/CosmicMutations_AtInformativeGenes_cosmic_v56_151111.3.bed'
genelist='/exds/users/oabalbin/projects/exomes/genelists/informative_genes.txt'
read_cosmic_mutation_file(ifile,ofile, list_of_names(open(genelist), pfield))
        
    
